% PRODUCEFIG3.M
%
% Produces figure 3 in the paper.

clear all;
%% Initial Parameters
% Initial Structure (holds up to period ta+T) 
% Monetary policy
target_inflation = 5;                     % Target rate of annual inflation (per cent);
rho_r            = 0.65;                  % Persistence of the interest rate in the Taylor-rule.
rho_pi           = 0.5;                   % Reaction to deviations of inflation from pitarget.
rho_y            = 0.1;                   % Reaction to deviation of output from steady state. 
rho_g            = 0.2;                   % Reaction to output growth
% Households and firm's 
z     = 1.1;                                  % Steady state level of TFP.
beta  = 0.9925;                               % Discount factor.
sigma = 1;                                    % equals to 1 for log-utility on consumption; 
alpha = 0.25;                                 % Share of `backward-looking' in pricing.
psi   = 0.1;                                  % Captures the degree of nominal rigidities: psi = (theta - 1)/ phi. 
% Exogenous process
rho_a = 0.9;                                  % Persistence of demand shock
rho_e = 0.9;                                  % Persistence of mark-up shock
rho_z = 0.9;                                  % Persistence of technology shock
siga  = 0.020;                                % Std deviation of innovation to the demand process.
sige  = 0.001;                                % Std deviation of innovation to the mark_up process. 
sigz  = 0.007;                                % Std deviation of innovation to the technology process. 
sigr  = 0.002;                                % Std deviation of innovation to the interest rate. 

Omega = [siga  0   0    0;
          0  sige  0    0;
          0    0  sigz  0;
          0    0    0   sigr];
% Other stats
pitarget = 0.0125; % target_inflation = 400*pitarget  
imp = Omega;


up_to = 16;    % Length of the simulation
ta = 4;        % Time of the announcement
T = 4;         % Implementation of the change at T + ta

time = 1:up_to;
% Parameter changes
% Monetary policy
Dtarget_inflation = -2.5; % In annual percent inflation
Drho_r            = 0; 
Drho_pi           = 0;
Drho_y            = 0;
Drho_g            = 0;
% Households and firm's 
Dz                = 0;
Dbeta             = 0;
Dsigma            = 0;
Dalpha            = 0;
Dpsi              = 0;
% Exogenous process
Drho_a            = 0;
Drho_e            = 0;
Drho_z            = 0;
Dsiga             = 0;
Dsige             = 0;
Dsigz             = 0;
Dsigr             = 0;

% Final Parameters
% New Structure (holds from period ta+T+1 onwards) 
target_inflation_n = target_inflation + Dtarget_inflation  ; 
rho_r_n = rho_r + Drho_r; 
rho_pi_n = rho_pi + Drho_pi; 
rho_y_n = rho_y + Drho_y;  
rho_g_n = rho_g + Drho_g; 
z_n = z + Dz; 
beta_n = beta + Dbeta; 
sigma_n = sigma + Dsigma; 
alpha_n = alpha + Dalpha; 
psi_n = psi + Dpsi; 
rho_a_n = rho_a + Drho_a; 
rho_e_n = rho_e + Drho_e; 
rho_z_n = rho_z + Drho_z; 
siga_n =  siga + Dsiga;  
sige_n =  sige + Dsige; 
sigz_n =  sigz + Dsigz; 
sigr_n =  sigr + Dsigr; 

Omega_n = [siga_n  0   0    0;
          0  sige_n  0    0;
          0    0  sigz_n  0;
          0    0    0   sigr_n];

% Other stats
pitarget_n = pitarget/2;
rss_n = (1+pitarget_n)/beta_n; % Steady state nominal interest rate.
rss_annual_n = ((1+pitarget_n)/beta_n)^4-1; % Steady state nominal interest rate in annual terms. 


%% Response under Initial Structure 
% Sets up the Stationary version of Ireland's (2007) New Keynesian model that we use for the applications in our paper. 
%
%                        1       2     3     4         5              6          7       8        9       10 
% Variables in Y(t) = [ yhat(t) pi(t) r(t) ghat(t) E_t(yhat(t+1)) E_t(pi(t+1)) ahat(t) ehat(t) ln(z(t)) zhat(t)]'
% 
%
%                           1      2    3     4
% Variables in eps(t) = [ ea(t) ee(t) ez(t) er(t) ]'
%
%
% Variables in eta(t) = [etay(t) etapi(t)]'

% Load Dimensions
n1 = 6;
n2 = 2;
nleads = 1;
s = nleads;
nfcast = 2;
nshocks = 4;
nvars = n1+n2+nfcast;
n = nvars;
k = nfcast;

% Load Initial Matrices
% Note: row indicates equation number and column indicates variables number
GAM0 = zeros(nvars);
GAM0(1,1) = 1; 
GAM0(1,3) = 1/sigma;
GAM0(1,5) = -1;
GAM0(1,6) = -1/sigma; 
GAM0(1,7) = -(1-rho_a)/sigma;
GAM0(2,1) = -psi*sigma;
GAM0(2,2) = 1+beta*alpha;
GAM0(2,6) = -beta;
GAM0(2,8) = 1;
GAM0(2,10) = psi;
GAM0(3,1) = -rho_y;
GAM0(3,2) = -rho_pi;
GAM0(3,3) = 1;
GAM0(3,4) = -rho_g;
GAM0(4,7) = 1;
GAM0(5,8) = 1;
GAM0(6,9) = 1;
GAM0(7,1) = -1;
GAM0(7,4) = 1;
GAM0(8,9) = -1;
GAM0(8,10) = 1;
GAM0(9,1) = 1;
GAM0(10,2) = 1;

GAM1 = zeros(nvars);
GAM1(2,2) = alpha;
GAM1(3,3) = rho_r;
GAM1(4,7) = rho_a;
GAM1(5,8) = rho_e;
GAM1(6,9) = rho_z;
GAM1(7,1) = -1;
GAM1(9,5) = 1;
GAM1(10,6) = 1;

C = zeros(nvars,1);
C(1,1) = -(1/sigma)*log(beta);
C(2,1) = (1+beta*alpha-alpha-beta)*pitarget;
C(3,1) = (1-rho_r)*(pitarget-log(beta)) - rho_pi*pitarget;
C(6,1) = (1-rho_z)*log(z);
C(8,1) = -log(z);

PSI = zeros(nvars,nshocks);
PSI(3,4) = 1;
PSI(4,1) = 1;
PSI(5,2) = 1;
PSI(6,3) = 1;

PPI = zeros(nvars,nfcast);
PPI(9,1) = 1;
PPI(10,2) = 1;

[S0_i, S1_i, S2_i, S3_i, unique_i] = Smats(C, GAM0, GAM1, PSI, PPI);
y0 = inv(eye(n) - S1_i)*S0_i;

yirfini = zeros(n,up_to);

% Impulse response to ea(t)
yirfini(:,1) = inv(eye(n)-S1_i)*S0_i + S2_i*imp(:,1);
for t = 1:(up_to-1)
    yirfini(:,t+1) = S0_i + S1_i*yirfini(:,t);
end

%% Response under Final Structure 
GAM0n = zeros(nvars);
GAM0n(1,1) = 1; 
GAM0n(1,3) = 1/sigma_n;
GAM0n(1,5) = -1;
GAM0n(1,6) = -1/sigma_n; 
GAM0n(1,7) = -(1-rho_a_n)/sigma_n;
GAM0n(2,1) = -psi_n*sigma_n;
GAM0n(2,2) = 1+beta_n*alpha_n;
GAM0n(2,6) = -beta_n;
GAM0n(2,8) = 1;
GAM0n(2,10) = psi_n;
GAM0n(3,1) = -rho_y_n;
GAM0n(3,2) = -rho_pi_n;
GAM0n(3,3) = 1;
GAM0n(3,4) = -rho_g_n;
GAM0n(4,7) = 1;
GAM0n(5,8) = 1;
GAM0n(6,9) = 1;
GAM0n(7,1) = -1;
GAM0n(7,4) = 1;
GAM0n(8,9) = -1;
GAM0n(8,10) = 1;
GAM0n(9,1) = 1;
GAM0n(10,2) = 1;

GAM1n = zeros(nvars);
GAM1n(2,2) = alpha_n;
GAM1n(3,3) = rho_r_n;
GAM1n(4,7) = rho_a_n;
GAM1n(5,8) = rho_e_n;
GAM1n(6,9) = rho_z_n;
GAM1n(7,1) = -1;
GAM1n(9,5) = 1;
GAM1n(10,6) = 1;

Cn = zeros(nvars,1);
Cn(1,1) = -(1/sigma_n)*log(beta_n);
Cn(2,1) = (1+beta_n*alpha_n-alpha_n-beta_n)*pitarget_n;
Cn(3,1) = (1-rho_r_n)*(pitarget_n-log(beta_n)) - rho_pi_n*pitarget_n;
Cn(6,1) = (1-rho_z_n)*log(z_n);
Cn(8,1) = -log(z_n);

PSIn = zeros(nvars,nshocks);
PSIn(3,4) = 1;
PSIn(4,1) = 1;
PSIn(5,2) = 1;
PSIn(6,3) = 1;

PPIn = zeros(nvars,nfcast);
PPIn(9,1) = 1;
PPIn(10,2) = 1;

[n l] = size(PSI);

% Reduced form solution matrices of final structure.
[S0_f, S1_f, S2_f, S3_f, unique_f] = Smats(Cn, GAM0n, GAM1n, PSIn, PPIn);
y0_new = inv(eye(n) - S1_f)*S0_f;

% Impulse response to ea(t)
yirffin(:,1) = inv(eye(n)-S1_f)*S0_f + S2_f*imp(:,1);
for t = 1:(up_to-1)
    yirffin(:,t+1) = S0_f + S1_f*yirffin(:,t);
end

%% Response with T = 4
% Compute transition period
% Store endogenous variables in the transition period
y_transition_4 = zeros(n,T);
% Load unanticipated shocks (use in stochastic simulations)
et = zeros(l,2*T); % size
% Store variables for overall response
yirfT4 = zeros(n,up_to);

% Construct sequences of time-varying matrices 
GAM0s = zeros(n*(T+1),n);
GAM1s = zeros(n*(T+1),n);
Cs = zeros(n*(T+1),1);
PSIs = zeros(n*(T+1),l);
PPIs = zeros(n*(T+1),k);

GAM0s(1:n*T,:) = repmat(GAM0,T,1);
GAM0s(n*T+1:end,:) = GAM0n;

GAM1s(1:n*T,:) = repmat(GAM1,T,1);
GAM1s(n*T+1:end,:) = GAM1n;

Cs(1:n*T,:) = repmat(C,T,1);
Cs(n*T+1:end,:) = Cn;

PSIs(1:n*T,:) = repmat(PSI,T,1);
PSIs(n*T+1:end,:) = PSIn;

PPIs(1:n*T,:) = repmat(PPI,T,1);
PPIs(n*T+1:end,:) = PPI;

complex_tol = 1e-6;
eps_a = zeros(l,T); % Load anticipated shocks

% IRF up to period ta-1
yirfT4(:,1:ta-1) = yirfini(:,1:ta-1);

% For initial period of the transition
[yseq, S0, S1, S2, Sy, Sf, Sz, A, B] = LREAnt(n1,n2,s,T,GAM0s,GAM1s,Cs,PSIs,PPIs,yirfT4(:,ta-1),eps_a, et(:,ta),complex_tol);
y0 = yseq(:,1);
y_transition_4(:,1) = yseq(:,1);
% Loop is useful for stochastic simulation.
for t = 2:T
    eps_a = zeros(l,T-t+1);
    GAM0s(1:n,:) = [];
    GAM1s(1:n,:) = [];
    Cs(1:n,:) = [];
    PSIs(1:n,:) = [];
    PPIs(1:n,:) = [];
   [yseq, S0, S1, S2, Sy, Sf, Sz, A, B] = LREAnt(n1,n2,s,T-t+1,GAM0s,GAM1s,Cs,PSIs,PPIs,y0,eps_a, et(:,ta+t-1),complex_tol);
    y0 = yseq(:,1);
    y_transition_4(:,t) = yseq(:,1);
end

yirfT4(:,ta:T+ta-1) = y_transition_4;

% End of Response with final structure

for t = T+ta:up_to
     yirfT4(:,t) = S0_f + S1_f * yirfT4(:,t-1);
end


%% Transformation of Units
output_resp = yirfT4(1,:);
g_resp = 400*yirfT4(4,:);
inflation_resp = yirfT4(2,:)*400;
r_resp = yirfT4(3,:)*400;

output_resp_ini = yirfini(1,:);
g_resp_ini = 400*yirfini(4,:);
inflation_resp_ini = yirfini(2,:)*400;
r_resp_ini = yirfini(3,:)*400;

output_resp_fin = yirffin(1,:);
g_resp_fin = 400*yirffin(4,:);
inflation_resp_fin = yirffin(2,:)*400;
r_resp_fin = yirffin(3,:)*400;

%% Plot Options
lw = 2;    % LineWidth
figure('Name','Figure 3 in the paper: Demand shock -- ea(t) shock 1','NumberTitle','off')
subplot(2,2,1), plot(time,output_resp_ini,'b--',time,output_resp_fin,'r--',time,output_resp,'g','LineWidth',lw), 
                title('Output','fontsize',14), ylabel('%'), legend('Initial Structure','Final Structure','Actual Path'),grid,set(gca,'xtick',[ta, T+ta])
subplot(2,2,2), plot(time,inflation_resp_ini,'b--',time,inflation_resp_fin,'r--',time,inflation_resp,'g','LineWidth',lw),
                title('Inflation','fontsize',14), ylabel('%'),legend('Initial','Final','Path'),grid,set(gca,'xtick',[ta, T+ta])
subplot(2,2,3), plot(time,r_resp_ini,'b--',time,r_resp_fin,'r--',time,r_resp,'g','LineWidth',lw),
                title('Nominal Interest Rate','fontsize',14), ylabel('%'),legend('Initial','Final','Path'),grid,set(gca,'xtick',[ta, T+ta])  
subplot(2,2,4), plot(time,g_resp_ini,'b--',time,g_resp_fin,'r--',time,g_resp,'g','LineWidth',lw),
                title('Output Growth','fontsize',14), ylabel('%'),legend('Initial','Final','Path'),grid,set(gca,'xtick',[ta, T+ta])
orient landscape


